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Pore  network  simulations  are  performed  to  study  water  transport  in  a  model  gas  diffusion  layer  (GDL) 
of  polymer  electrolyte  membrane  fuel  cells  (PEMFCs)  in  relation  with  the  change  in  hydrophobicity  that 
might  be  due  to  aging  or  temperature  effect.  The  change  in  hydrophobicity  is  taken  into  account  by 
changing  randomly  the  fraction  of  hydrophilic  elements,  pores  or  throats,  in  the  network.  The  transport 
and  equilibrium  properties  of  the  model  GDL  are  computed  as  a  function  of  liquid  saturation  as  well  as 
at  breakthrough  varying  the  fraction  of  hydrophilic  elements.  The  results  indicate  that  the  hydrophilic 
element  percolation  threshold  marks  the  transition  between  two  domains.  The  system  is  found  to  be 
weakly  dependent  on  the  fraction  of  hydrophilic  elements  as  long  as  this  fraction  is  below  the  percolation 
threshold  whereas  an  increase  in  wettability  above  the  percolation  threshold  favours  a  greater  blockage 
of  the  pore  space  by  the  water  and  therefore  a  diminished  access  of  gas  to  the  catalyst  layer.  This  model 
may  help  assess  the  effect  of  a  change  in  wettability  on  the  fuel  cell  performance  and  may  also  help 
suggest  better  GDL  designs  in  relation  with  the  water  management  problem  in  PEMFCs. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  continues  to  be  a  major  issue  for  proton 
exchange  membrane  fuel  cells,  e.g.  [1]  and  references  therein,  and 
it  is  now  clear  that  the  gas  diffusion  layer  (GDL)  is  an  important 
element  in  the  prospect  of  a  better  water  management.  Optimiz¬ 
ing  a  GDL  for  the  water  management  problem  is,  loosely,  having 
enough  water  to  ensure  the  good  hydration  of  the  electrolyte  while 
minimizing  the  impact  of  liquid  water  in  the  pores  on  the  gas 
access  to  the  catalyst  layer.  Despite  numerous  studies  devoted  to 
GDL,  e.g.  [2]  and  references  therein,  there  is  still  no  obvious  way 
of  optimizing  the  fuel  cell  performances  by  playing  with  the  GDL 
properties.  One  reason  is  that  the  effect  of  GDL  wettability  prop¬ 
erties  is  still  not  fully  well  understood.  Pore  network  simulations, 
e.g.  [3],  suggest  that  a  sufficiently  hydrophobic  GDL  is  better  than 
a  hydrophilic  one.  This  is  because  water  invasion  in  a  hydrophobic 
porous  medium  in  the  quasi-static  limit  that  is  expected  to  pre¬ 
vail  for  flows  in  fuel  cells  leads  to  capillary  fingering  whereas  the 
invasion  pattern  is  compact  in  a  hydrophilic  medium.  The  capillary 
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fingering  pattern  leaves  many  pores  not  invaded  by  the  liquid  and 
therefore  available  for  the  gas  transport.  In  contrast,  the  compact 
invasion  pattern  blocks  rapidly  the  gas  transport  across  the  GDL. 
This  sounds  consistent  with  the  fact  that  a  GDL  is  generally  treated 
with  a  hydrophobic  fluoropolymer.  However,  the  distribution  of  the 
fluoropolymer  on  the  carbon  fibres  forming  the  GDL  microstruc¬ 
ture  is  heterogeneous,  e.g.  [1]  and  references  therein.  The  SEM 
micrographs  of  GDL  reported  in  [  1  ]  indicate  that  there  are  more  flu¬ 
oropolymer  at  the  junctions  between  fibres  for  example  and  that 
other  parts  of  the  fibres  might  be  not  coated  by  the  fluoropoly¬ 
mer.  However  the  more  precise  ToF-SIMS  images  also  reported  in 
[1  ]  lead  to  the  conclusion  that  the  fluoropolymer  is  in  fact  present 
everywhere  on  the  pore  wall.  Hence  Teflon  would  be  present  in 
the  form  of  thin  films  even  in  the  regions  of  the  microstructure 
which  appear  as  seemingly  not  teflonized  in  the  SEM  micrographs. 
Since  wetting  is  generally  dominated  by  the  top  few  monolayers 
of  a  surface  [4],  it  is  anticipated  that  the  wetting  properties  of  a 
teflonized  GDL  microstructure  are  uniform  and  close  to  the  ones  of 
water  on  Teflon.  This  leaves  us  with  two  options.  The  first  option 
is  to  consider,  as  suggested  in  [1]  as  well  as  in  [5],  that  the  teflo- 
nization  process  is  effective  and  therefore  that  the  wettability  of 
a  GDL  is  spatially  uniform,  that  is  the  (advancing)  contact  angle 
varies  locally  within  a  quite  narrow  range  around  115°  or  possi¬ 
bly  more.  The  second  option,  considered  for  example  in  [6-9],  is  to 
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consider  that  a  GDL  should  be  envisioned  as  a  system  of  mixed  wet¬ 
tability,  that  is  some  regions  of  the  microstructure  are  hydrophilic 
(with  the  local  contact  angle  0  expected  to  be  on  the  order  of  that 
of  water  on  carbon,  0~8O°)  whereas  others  are  hydrophobic,  i.e. 
coated  with  the  fluoropolymer.  How  to  distribute  the  hydrophilic 
and  hydrophobic  regions  on  the  pore  wall  surface  is  not  obvious 
since  we  do  not  know  for  example  the  size  of  the  hydrophilic  (resp. 
hydrophobic)  regions.  As  for  example  in  [6],  the  simple  option  con¬ 
sidered  in  this  paper  is  to  suppose  that  a  fraction /of  the  pores  is 
hydrophilic,  the  others  (corresponding  to  the  fraction  1  -f)  being 
hydrophobic  and  to  explore  the  influence  of/on  water  invasion  or 
water  distribution  in  the  network.  It  could  be  argued  that  it  is  some¬ 
what  not  consistent  to  consider  mixed  wettability  systems  since  the 
recent  results  presented  in  [1,5]  indicate  that  it  seems  more  appro¬ 
priate  to  assume  a  uniform  wettability.  It  is  shown,  however,  that 
our  mixed  wettability  model  is  in  good  agreement  with  the  most 
recent  measurements  of  the  capillary  pressure  curve.  Thus  it  can¬ 
not  be  readily  concluded  from  these  measurements  that  the  GDL 
is  of  uniform  wettability.  Even  if  the  assumption  of  spatially  uni¬ 
form  wettability  is  the  correct  assumption  for  the  current  teflonized 
GDLs,  it  remains  interesting,  anyway,  to  explore  whether  systems  of 
mixed  wettability  would  be  better  for  the  water  management  issue. 
Also,  it  is  widely  admitted  that  the  GDLs  lose  hydrophobicity  dur¬ 
ing  the  course  of  a  PEMFC  lifetime,  e.g.  [10].  It  is  unclear,  however, 
whether  this  change  in  hydrophobicity  occurs  everywhere  within 
the  pore  space  or  rather  according  to  a  scenario  where  the  change 
in  hydrophobicity  would  be  due  to  the  increase  in  the  fraction  of 
pore  wall  surface  becoming  hydrophilic  over  time.  Also,  it  has  been 
reported  that  a  GDL  tends  to  become  more  hydrophilic  after  expo¬ 
sure  to  liquid  water  at  temperatures  representative  of  the  PEMFC 
operating  conditions  [1 1,12].  Increasing  the  fraction  of  hydrophilic 
pores  in  a  mixed  wettability  model  is  a  convenient  way  to  explore 
the  consequences  of  this  change  in  wettability.  As  shown  in  [3],  see 
also  [13]  and  references  therein,  it  is  possible  to  explore  numeri¬ 
cally  the  effect  of  a  spatially  uniform  change  in  wettability  in  the 
range  of  intermediate  contact  angles  expected  for  the  GDLs,  i.e.  typ¬ 
ically  in  the  “critical”  range  [70-115°],  in  model  two-dimensional 
systems.  However,  a  three-dimensional  version  of  this  model  is 
not  available  and  direct  simulations  at  the  pore  scale  based  on  lat¬ 
tice  Boltzmann  methods,  e.g.  [14],  or  Navier-Stokes  equations  are 
too  much  computational  time  consuming  for  a  convenient  explo¬ 
ration  of  the  effect  of  the  contact  angle.  Thus,  the  most  convenient 
way  at  the  moment  for  exploring  a  change  in  wettability  within 
the  framework  of  three-dimensional  PNM  is  to  vary  the  fraction  of 
hydrophilic  pores. 

The  paper  is  organized  as  follows.  The  pore  network  model  is 
presented  in  Section  2.  In  Section  3,  we  discuss  in  more  depth  the 
wettability  of  GDL  from  comparisons  between  pore-network  sim¬ 
ulations  and  measured  capillary  pressure  curves.  The  computation 
of  the  effective  diffusion  coefficient  and  of  the  relative  perme¬ 
abilities,  two  important  parameters  for  the  continuum  models,  is 
presented  in  Section  4.  Then,  we  present  in  Section  5  simulations 
of  the  water  invasion  of  a  GDL  considering  various  wettability  sce¬ 
narios  so  as  to  explore  the  influence  of  wettability  on  pore  blockage 
by  the  water.  Note  that  the  structural  changes  (pore  sizes,  poros¬ 
ity,  etc.)  that  might  result  from  varying  the  load  in  PTFE  are  not 
considered. 


2.  Pore  network  model  (PNM) 

The  PNMs  are  based  on  the  representation  of  the  pore  space 
in  terms  of  a  network  of  pores  (or  sites)  connected  by  throats 
(or  bonds).  The  “pores”  roughly  correspond  to  the  larger  voids 
whereas  the  throats  connecting  the  pores  correspond  to  the  con¬ 
strictions  of  the  pore  space.  In  the  most  advanced  developments, 


see  for  example  [15]  and  references  therein,  the  pore  network  is 
constructed  from  direct  imaging,  usually  by  micro-X-ray  comput¬ 
erized  tomography.  An  alternative  is  to  use  synthetic  3D  structure 
generated  numerically.  This  leads  to  “morphological”  pore  net¬ 
works  since  the  pore  network  is  constructed  directly  from  the 
“real”  microstructure.  The  method  is  for  example  illustrated  in 
[3]  with  a  two-dimensional  network  constructed  from  a  model 
fibrous  medium.  Although  using  morphological  pore  networks  is 
certainly  the  most  attractive  approach,  using  simpler  pore  net¬ 
works  can  be  very  instructive.  In  fact,  morphological  pore  networks 
have  rarely  been  used  in  relation  with  fuel  cell  related  problems, 
see  however  [16].  Most  of  the  studies,  e.g.  [3,17-21],  are  based  on 
simpler  networks,  such  as  cubic  networks  in  3D  or  square  networks 
in  2D.  Computations  are  easier  than  for  a  morphological  network 
and  it  is  still  possible  to  incorporate  information  from  the  “real” 
microstructure  such  as  throat  size  and  pore  size  distributions.  Also, 
to  understand  some  effects,  a  simple  network  is  generally  sufficient. 
Accordingly  with  the  above,  a  simple  cubic  network  is  considered 
in  what  follows.  As  discussed  in  several  places  in  the  paper  from 
comparisons  with  experimental  data,  this  simple  model  seems  to 
be  not  always  sufficient,  however,  to  represent  fibrous  materials 
accurately.  The  development  of  advanced  pore  network  models  for 
fibrous  materials  is  in  fact  an  open  area  of  research.  We  believe, 
however,  that  the  model  considered  in  the  present  study  is  suf¬ 
ficient  to  explore  the  effect  of  a  mixed  wettability,  which  is  the 
primary  objective  of  the  study. 

As  sketched  in  Fig.  1,  pores  of  cubic  shape  are  regularly  placed 
on  a  3D  cartesian  grid  (with  a  denoting  the  lattice  spacing,  that 
is  the  distance  between  two  adjacent  pores).  Two  first  neigh¬ 
bour  pores  are  linked  by  a  channel  of  square  cross-section.  Such 
a  channel  is  referred  to  as  a  bond  or  throat.  The  pore  size  dp , 
which  corresponds  to  the  diameter  of  the  largest  sphere  inscribed 
within  the  pore,  is  randomly  distributed  according  to  a  probability 
law  in  the  range  [dpm[n,  dp  max].  As  in  [17],  a  Weibull  distribu¬ 
tion  is  used  in  the  present  study.  More  precisely,  the  pore  sizes 
are  randomly  specified  using  the  expression  dp  =  dpm[n  +  (dpmax  - 


up  min ) 


i  i/y 


with  5  =  0.1, 


a  [{-5  ln(X(l  -  exp(-l ./ <5))  +  exp(-l  ./(5))}1 

y  =  4.7;  X  is  the  random  number  drawn  in  the  interval  [0,1].  The 
size  dt  of  a  throat  is  then  specified  as  equal  to  the  minimum  size 
of  its  two  adjacent  pores,  dt  =  min  (dpi,dpj).  Simulations  presented  in 
what  follows  have  been  performed  with  a  =  25.2  p,m,  dpmax  =  25  p,m, 
dpmin  =  10  |jim.  These  values  are  representative  of  pore  size  distri¬ 
butions  in  GDL,  e.g.  [17],  and  lead  to  an  average  porosity  s  of  0.77 
for  our  model  porous  medium. 

The  network  size  is  characterized  by  the  number  of  pores  placed 
in  each  direction.  We  consider  a  40  x  40  x  1 0  network.  As  discussed 
in  [21],  this  size  is  considered  as  representative  of  a  GDL  unit  cell 
with  the  in-plane  extension  (~1  mm)  corresponding  roughly  to  the 
distance  between  two  channels  of  the  bipolar  plate.  The  pore  net¬ 
work  is  fully  hydrophobic  (resp.  hydrophilic)  when  all  throats  and 
pores  are  hydrophobic  (resp.  hydrophilic).  In  a  mixed  wettability 
pore  network  a  fraction  /  of  the  pores  and  of  the  throat  are  ran¬ 
domly  specified  as  hydrophilic  (0~8O°).  The  remaining  pores  and 
throats  are  hydrophobic  (0~  115°).  There  are  therefore  a  fraction 
1  -/of  hydrophobic  pores  and  throats  in  the  network.  For  conve¬ 
nience,  such  a  network  will  be  referred  to  as  a  network  containing  a 
fraction/of  hydrophilic  “pores”  but  this  means  in  fact  that  the  frac¬ 
tion  of  hydrophilic  throats  and  the  fraction  of  hydrophilic  pores  are 
both  equal  to/.  No  correlation  is  introduced  between  the  wettabil¬ 
ity  of  a  pore  and  that  of  adjacent  throats.  A  realization  of  network 
is  obtained  by  assigning  pore  and  throat  size  randomly  and  0  (80° 
or  110°)  randomly  for  a  given  hydrophilic  fraction.  All  the  results 
presented  in  the  following  are  averaged  over  100  realizations  (the 
results  presented  in  [22]  indicate  that  this  is  a  sufficient  number  of 
realizations  to  obtain  the  average  behaviours). 
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Fig.  1.  Schematic  of  cubic  pore  network.  The  size  of  network  is  40  x  40  x  10.  As  shown  in  the  figure,  the  GDL  thickness  is  252  (Jim  and  there  are  10  pores  25.2  p,m  apart  across 
the  GDL. 


3.  Local  invasion  rules 


It  is  widely  admitted  that  liquid  invasion  in  a  GDL  can  be  com¬ 
puted  in  the  quasi-static  limit,  that  is  the  invasion  is  controlled 
by  capillary  effects,  see  [3]  or  [22]  for  more  details.  When  the 
computation  of  properties  is  performed  in  relation  with  contin¬ 
uum  models,  the  classical  assumption  of  local  capillary  equilibrium 
is  assumed  and  this  corresponds  again  to  quasi-static  water  dis¬ 
tribution  in  the  pore  space.  Accordingly,  the  influence  of  viscous 
effect  is  ignored  in  this  paper  and  only  capillary  effects  are  taken 
into  account.  When  the  porous  matrix  is  sufficiently  hydropho¬ 
bic,  quasi-static  water  invasion  can  be  simulated  rather  simply 
using  the  well-known  invasion  percolation  (IP)  algorithm  [23]. 
When  the  medium  is  hydrophilic,  the  computation  is  significantly 
more  involved  because  of  cooperative  mechanisms  controlling  the 
growth  of  the  interface  between  the  two  fluids  within  the  pore 
space  [3,21  ],  and  references  therein.  A  relatively  simple  but  approx¬ 
imate  algorithm  was,  however,  proposed  in  [24],  see  also  [25].  A 
similar  algorithm  is  used  in  the  present  study.  Note  that  this  algo¬ 
rithm  becomes  identical  to  the  IP  algorithm  when  the  network  is 
fully  hydrophobic  (f=  0).  Hence  all  the  simulations  are  performed 
using  the  local  expressions  given  in  this  section.  The  capillary  pres¬ 
sure  is  conventionally  defined  as 


Pc  =Pw-  Pg 


(1) 


i.e.  as  the  difference  between  the  pressure  in  the  water  and  the 
pressure  in  the  gas  phase.  Then  we  define  the  capillary  pressure 
threshold  pcth  associated  with  each  element,  pore  or  throat,  of 
the  pore  network.  The  capillary  pressure  threshold  of  a  throat  is 
expressed  using  the  Young-Laplace  equation  as 


Pcth  =  - 


2(7  cos  0 
rt 


+  Pi 


(2) 


where  rt  &  dt\ 2;  Pi  is  a  numerical  constant  with  P\  =  0  for  0  =  1 1 5° 
and  P\  =5500  for  0  =  80°.  The  motivation  for  the  introduction  of 
constant  Pi  is  given  in  the  next  section  (Section  4).  Eq.  (2)  is  only 
an  approximation.  One  can  refer  for  example  to  [26]  for  a  more 
accurate  evaluation  of  the  capillary  pressure  threshold  of  a  throat 
of  square  cross-section.  Eq.  (2)  leads  to  the  correct  hierarchy  in 
terms  of  capillary  pressure  threshold  and  this  is  sufficient  for  the 


pore  network  simulations  (PNS).  The  contact  angle  0  is  measured 
in  the  denser  fluid  (water);  a  is  the  surface  tension.  The  capillary 
pressure  threshold  of  a  pore  depends  on  the  number  n  of  adjacent 
throats  filled  with  the  less  dense  fluid  (gas  phase)  and  reads 


Pcthp  = 


2(7  COS0 

rn 


+  Pi 


where 


n 

Y)b<x>rt> 


i=i 


(3) 


(4) 


where  dp  is  the  pore  size,  rti  is  the  half-size  of  a  adjacent  throat  filled 
with  air,  x,-  is  a  random  number  drawn  in  the  range  0-1 .  The  coeffi¬ 
cients  b\  are  defined  as:  b\  =0.,  b2  =  5 ./ft,  b-s,  =  10. /ft,  b4  =  50. /ft, 
b5  =  100./rt;  ft  is  the  mean  throat  radius  (mean  of  the  pore  sur¬ 
rounding  throat  radii).  When  n  is  small,  many  surrounding  throats 
are  already  filled  with  water  and  water  invasion  of  the  pore  is 
favoured  (when  0<9O°).  As  discussed  in  [27],  Eqs.  (2)-(4)  define 
invasion  potentials.  The  element,  pore  or  throat,  of  lowest  poten¬ 
tial  (=lowest  capillary  pressure  threshold)  is  invaded  at  each  step 
of  invasion  to  simulate  a  quasi-static  invasion  (see  Section  7).  The 
coefficients  b{  are  specified  so  as  to  obtain  a  compact  pattern  in 
a  100%  hydrophilic  pore  network.  The  values  of  these  coefficients 
indicated  above  lead  to  the  desired  pattern.  One  can  refer  to  the 
Appendix  in  [27]  for  more  details  on  the  invasion  potential  hier¬ 
archy  needed  to  simulate  adequately  the  quasi-static  invasion  of  a 
porous  medium  by  a  wetting  fluid. 

As  explained  in  the  next  section,  the  constant  Pi  is  added  to  Eqs. 
(2)  and  (3)  to  generate  positive  capillary  pressures  in  fully  or  par¬ 
tially  hydrophilic  networks  so  as  to  be  in  agreement  with  available 
experimental  data.  This  might  appear  as  quite  artificial.  However,  it 
should  be  realized  again  that  the  key  point  is  that  Eqs.  (2)-(4)  lead 
to  a  correct  hierarchy  in  terms  of  local  invasion  capillary  thresh¬ 
olds,  i.e.  for  example  a  hydrophilic  element  is  always  preferentially 
invaded  compared  to  a  hydrophobic  one.  Adding  the  constant  Pi 
does  not  affect  this  hierarchy  and  therefore  the  simulated  water 
distributions  and  consequently  the  results  presented  in  what  fol¬ 
lows,  except  as  regards  the  capillary  pressure  curve  for  which  a 
much  better  agreement  is  obtained  by  introducing  the  constant  Pi . 


1148 


S.P.  Kuttanikkad  et  al.  /  Journal  of  Power  Sources  196(2011)  1 145-1155 


f=0  (all  pores  hydrophobic)  f =0.2  <f.  f =0.5  <f.  f =0.6  <fc 


f=0.7~fc  f=0.8>fc  f=0.9  >fc  /=/(all  pores  hydrophilic) 


Fig.  2.  Invasion  patterns  for  various  values  of  the  hydrophilic  pore-throat  fraction/ in  a  2D  50  x  50  square  network.  The  invading  fluid  (water  in  white)  is  injected  from  the 
bottom  using  the  algorithm  described  in  Section  7.  The  defending  fluid  (in  black)  escapes  from  the  top.  A  capillary  fingering  pattern  is  observed  for/~  0.  A  compact  pattern 
is  obtained  for/~l.  The  pattern  obtained  for  /~/c  resembles  the  capillary  fingering  pattern  but  actually  results  from  the  fact  the  percolating  subnetwork  formed  by  the 
hydrophilic  pores  and  throats  is  a  percolation  cluster.  An  invasion  percolation  cluster  0)  is  about  the  same  fractal  object  as  the  hydrophilic  pores  percolation  cluster 
This  is  why  both  cases  (f~  0  and /~/c)  lead  to  the  same  type  of  phase  distribution  despite  the  change  in  the  local  mechanisms  controlling  the  growth  of  the  liquid/gas 
interface  in  the  hydrophobic  regions  and  the  hydrophilic  ones  respectively.  For  this  bond-site  percolation  network,  fc  is  approximately  equal  to  0.7. 


The  two-dimensional  invasion  patterns  obtained  at  break¬ 
through  using  these  local  rules  in  a  square  network,  see  Section 
7  for  the  details  on  the  invasion  algorithm,  and  shown  in  Fig.  2, 
exemplify  the  effect  of  the  change  in  wettability.  As  can  be  seen, 
the  patterns  are  roughly  similar  (i.e.  ramified  and  not  compact) 
as  long  as  /  is  lower  than  the  so-called  percolation  threshold  fc. 
As  discussed  in  more  details  in  Section  7,  fc  is  the  critical  value 
above  which  at  least  one  continuous  path  of  hydrophilic  pores 
and  throats  exists  between  the  inlet  and  the  outlet.  The  pat¬ 
tern  becomes  more  and  more  compact  as  /  increases  above  fc.  In 
a  hydrophilic  network  (f=  1 ),  the  pattern  is  characterized  by  an 
almost  flat  travelling  front.  We  chose  to  present  two-dimensional 
patterns  in  Fig.  2  only  for  clarity  since  a  2D  pattern  is  much  easier 
to  represent  than  a  3D  pattern.  Analogous  patterns  are  obtained 
in  3D  as  well  and  this  should  be  kept  in  mind  when  discussing 
the  results  presented  in  the  next  sections.  All  the  results  pre¬ 
sented  in  the  remaining  of  the  paper  are  for  a  3D  40  x  40  x  10 
network. 

4.  Capillary  pressure  curve 

The  capillary  pressure  curve  Pc(5)  (or  equivalently  Pc(Sg )  where 
Sg  =  l  -S,  where  S  is  the  saturation  of  the  denser  phase  (=liquid 
water)  and  Sg  is  the  saturation  of  the  gas  phase)  is  one  of  the  key 
input  parameters  for  the  traditional  continuum  models.  Also,  the 
capillary  pressure  curve  gives  information  about  the  overall  wet¬ 
tability  of  a  porous  system,  e.g.  [1,5].  As  exemplified  for  example 
in  [6],  pore  network  models  can  be  used  to  determine  PC(S).  The 
capillary  curve  is  computed  using  an  algorithm  similar  to  the  one 
presented  in  [22],  but  which  now  takes  into  account  the  presence  of 
hydrophilic  pores.  This  algorithm  is  supposed  to  mimic  the  exper¬ 
imental  procedure  used  for  example  in  [1].  Due  to  the  random 
fluctuations,  PC(S)  has  to  be  determined  over  many  realizations 
of  the  network.  For  a  given  realization,  we  start  with  a  dry  net¬ 
work  and  determine  the  overall  capillary  pressure  as  a  function 
of  S  for  successive  states  of  hydrostatic  equilibrium  corresponding 
to  small  increment  dPw  in  the  liquid  water  pressure  (the  gaseous 


phase  pressure  is  kept  constant).  The  algorithm  used  to  determine 
the  saturation  evolution  right  after  a  pressure  increment  reads: 

1.  Identification  of  each  element  (pore  or  throat)  that  can  be 
invaded,  that  is  each  element  occupied  by  the  gaseous  phase 
such  that  its  capillary  pressure  threshold  pcth  is  lower  than  the 
imposed  pressure  difference  Pw  -  Pg  between  the  two  fluids. 

2.  Identification  of  clusters  formed  by  the  elements  identified  in  ( 1 ) 
using  a  cluster  identification  algorithm. 

3.  Identification  of  clusters  among  the  clusters  identified  in  (2)  in 
contact  with  the  invading  phase. 

4.  Invasion  of  all  clusters  identified  in  (3). 

5.  Computation  of  overall  saturation. 

Using  this  algorithm  leads  to  the  results  shown  in  Fig.  3.  The 
results  are  compared  with  the  experimental  results  presented  in 
[1]  in  Fig.  3a.  As  can  be  seen,  there  is  a  good  agreement  between 
the  simulations  and  the  experimental  results  over  a  wide  range  of 
intermediate  saturations  for  both  the  fully  hydrophobic  and  the 
fully  hydrophilic  network.  There  are  clearly  discrepancies  for  the 
low  and  high  saturations.  We  have  changed  the  pore  size  distri¬ 
bution  (PSD)  using  a  uniform  distribution  (each  pore  size  has  the 
same  probability)  and  observed  that  a  slightly  better  agreement 
could  be  obtained.  Thus  this  suggests  that  these  differences  could 
be  explained  in  part  by  the  choice  of  the  PSD.  For  example,  the  cap¬ 
illary  pressure  level  reached  in  the  experiments  (~14kPa)  is  much 
higher  than  in  the  simulations.  This  is  probably  due  to  the  fact  that 
very  small  pores  are  not  considered  in  the  simulations.  If  we  make 
the  fair  assumption  that  the  experimental  results  are  fully  reliable, 
the  comparison  shown  in  Fig.  3a  could  be  an  indication  that  the 
problem  comes  also  from  the  pore  network  model.  A  simple  cubic 
pore  network  as  considered  here  is  perhaps  not  sufficient  to  simu¬ 
late  accurately  the  capillary  pressure  curve  and  more  refined  pore 
network  models,  i.e.  typically  morphological  or  topological  pore 
network  models  as  mentioned  in  the  introduction,  would  deserve 
to  be  developed.  As  mentioned  before,  this  is  however  an  open 
area  of  research  as  regards  fibrous  materials  and  therefore  this  is 
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Fig.  3.  Computed  water  invasion  capillary  pressure  curve  for  various  fractions  of 
hydrophilic  pores  (solid  line  for  a  hydrophilic  pore  fraction  below  the  percolation 
threshold  fc,  dashed  line  above  fc ):  (a)  comparison  with  experimental  data:  the 
empty  circles  and  squares  correspond  to  the  measured  capillary  pressure  curves  for 
PTFE  treated  Toray  TGPH-090  (empty  circles)  and  untreated  Toray  TGPH-090  (empty 
squares)  (imported  and  adapted  from  [1  ]).  The  computations  were  performed  with 
a  contact  angle  of  1 1 5°  in  the  hydrophobic  pores  and  of  80°  in  the  hydrophilic  pores; 
(b)  influence  of  the  fraction  of  hydrophilic  pores/on  the  computed  capillary  pressure 
curve. 


left  for  a  future  work.  Here  we  consider  that  the  good  agreement 
obtained  over  a  wide  range  of  intermediate  saturations  is  sufficient 
to  explore  with  confidence  the  effect  of  a  mixed  wettability  with 
the  present  model. 

Another  aspect  of  the  experimental  results  deserves  to  be 
pointed  out.  Except  for  the  low  saturations  below  0.1,  the  exper¬ 
iments  lead  to  a  positive  capillary  pressure  (the  pressure  in  the 
liquid  water  is  greater  than  in  the  gas  phase)  for  both  the  treated 
and  untreated  GDLs.  This  was  expected  for  a  treated  GDL  since 
the  contact  angle  on  PTFE  is  significantly  greater  than  90°.  For  the 
untreated  GDL,  the  contact  angle  is  expected  to  be  lower  than  90° 
(0~8O°)  and  the  naive  consideration  of  the  fibrous  materials  as  a 
bundle  of  straight  tubes  leads  immediately  to  the  conclusion  that 
the  capillary  pressure  should  be  negative  (with  our  convention)  in 


a  hydrophilic  medium.  It  is  however  well  known  that  the  capillary 
pressure  threshold  of  a  constriction  depends  not  only  on  the  contact 
angle  but  also  on  the  constriction  geometry,  e.g.  [28].  The  consid¬ 
eration  of  a  more  realistic  geometry,  e.g.  [3,28],  leads  to  positive 
capillary  pressure  threshold.  To  mimic  this  effect,  we  have  chosen 
the  simpler  approach  consisting  to  add  a  positive  constant  Pi  in  Eqs. 
(2)  and  (3 )  for  the  contact  angle  0  =  80°.  As  can  be  seen  from  Fig.  3a,  a 
reasonable  agreement  with  the  experimental  data  is  obtained  with 
Pi  =  5500  in  Eqs.  (2)  and  (3).  The  model  being  adjusted  on  available 
capillary  pressure  experimental  data,  we  now  discuss  the  effect  of 
mixed  wettability  on  the  capillary  pressure  curve. 

As  can  be  seen  from  Fig.  3b,  the  results  for  a  fraction  of 
hydrophilic  pores  lower  than  about  50%  are  practically  identical  to 
those  for  a  fully  hydrophobic  network  except  for  the  high  satura¬ 
tions.  In  passing,  we  note  that  the  behaviour  for  the  high  saturations 
is  closer  to  the  experimental  results  (Fig.  3a  treated  GDL)  when 
/=  40%  for  example.  It  is  therefore  tempting  to  conclude  that  this 
is  an  indication  that  the  treated  GDL  is  in  fact  partially  hydropho¬ 
bic  but  again  more  advanced  pore  network  models  are  needed  to 
confirm  this  from  PNS.  Interestingly,  the  influence  of  Teflon  loading 
on  the  (water  invasion)  capillary  curve  is  found  to  be  quite  weak 
in  the  experiments  for  all  the  loading  tested  (up  to  40wt%),  e.g. 
[1].  A  similar  result  is  found  numerically  provided  that  the  frac¬ 
tion  of  hydrophilic  pores  is  lower  than  about  50%,  which  roughly 
corresponds  to  the  hydrophilic  pore/throat  network  percolation 
threshold  discussed  in  Section  7.  When  the  fraction  of  hydrophilic 
pores  is  greater  than  50%  the  results  are  markedly  different  as  can 
be  seen  from  Fig.  3. 

In  summary,  the  comparison  between  the  40  x  40  x  1 0  pore  net¬ 
work  simulations  and  the  experimental  results  presented  in  this 
section  shows  (i)  that  the  pore  network  model  can  simulate  rea¬ 
sonably  well  the  capillary  pressure  in  a  fully  hydrophobic  or  a  fully 
hydrophilic  system,  (ii)  that  the  assumption  of  mixed  wettability 
is  consistent  with  the  available  data  on  the  GDL  capillary  pressure 
curve  when  the  fraction  of  hydrophilic  pores  is  lower  than  about 
50%  and  (iii)  that  more  refined  pore  network  models  are  needed 
for  obtaining  a  still  better  match  with  the  experimental  data.  Our 
simulations  also  show  that  it  cannot  be  readily  concluded  from 
the  measurement  of  the  water  invasion  capillary  pressure  curve 
that  the  system  under  study  is  fully  hydrophobic  (as  suggested 
in  [1]  from  the  ToF-SIMS  data)  since  the  curve  is  not  sensitive, 
except  however  in  the  range  of  the  high  saturations,  to  a  change 
in  the  fraction  of  hydrophilic  pores  up  to  relatively  high  fractions 
of  hydrophilic  pores  (f<  50%). 


5.  Relative  permeabilities 

The  relative  permeabilities  are  also  important  parameters  of 
continuum  models.  In  this  section  we  discuss  the  influence  of  the 
fraction  of  hydrophilic  pores  on  the  relative  permeabilities.  The  pro¬ 
cedure  for  computing  the  intrinsic  permeability  and  the  relative 
permeabilities  is  described  in  several  papers,  e.g.  [22]  and  refer¬ 
ences  therein,  and  therefore  is  only  briefly  summarized  here.  Local 
hydraulic  conductances  are  assigned  to  each  throat  and  pore.  The 
same  expressions  as  the  ones  considered  in  [17]  are  used.  Express¬ 
ing  the  mass  conservation  equation  at  each  pore  yields  a  linear 
system  of  equations  that  is  solved  numerically  and  gives  the  pres¬ 
sure  at  each  node  (pores)  of  the  network.  As  boundary  conditions, 
arbitrary  pressures  Poutlet  and  Piniet are  prescribed  at  inlet  and  outlet 
surfaces  of  network  with  Pinlet  >  Poutiet-  This  yields  the  total  flow  Q. 
over  the  network.  Once  Q  is  determined,  the  permeability  I<  of  the 
network  is  found  from  Darcy’s  law.  For  determining  the  relative 
permeabilities  we  combine  the  procedure  described  in  Sec¬ 
tion  4  to  determine  the  fluids  distribution  within  the  network  and 
the  computation  of  permeability  described  above.  For  example  the 
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Fig.  4.  Through-plane  relative  permeability  of  liquid  water  (a,  b)  and  gas  (c,  d)  as  a 
hydrophilic  pore  fraction  below  the  percolation  threshold  fc,  dashed  line  above /c). 


effective  permeability  I<  I<rg  of  the  gas  phase  is  determined  using 
the  same  procedure  as  for  the  intrinsic  permeability  except  that 
only  the  sub-network  occupied  by  this  phase  is  considered.  Again 
more  details  on  this  type  of  computation  can  be  found  in  [17]. 

The  evolutions  of  I<r  for  our  network  are  shown  in  Fig.  4 
for  various  fractions  of  hydrophilic  pores  and  correspond  to  the 
through-plane  relative  permeability.  As  can  be  seen  from  Fig.  4,  the 
relative  permeability  curves  are  not  very  sensitive  to  the  fraction  of 
hydrophilic  pores  in  the  range  of/[0-50%].  For  a  given  intermediate 
value  of  the  saturation  (S  ~  0.5  for  example),  the  gaseous  phase  rela¬ 
tive  permeability  increases  slightly  with  the  fraction  of  hydrophilic 
pores  whereas  the  opposite  behaviour  is  observed  for  the  liquid 
water  relative  permeability.  As  can  be  seen  from  Fig.  4b  and  d,  the 
sensitivity  of  the  relative  permeabilities  to  the  hydrophilic  pore 
fraction  is  much  more  marked  for/>50%.  For  again  an  intermedi¬ 
ate  value  of  the  saturation,  5  ~  0.4-0.5,  the  gas  (resp.  liquid)  relative 
permeability  is  significantly  greater  (resp.  lower)  in  the  hydrophilic 
network  (f=  1 00%  in  Fig.  4).  As  discussed  in  the  next  Section,  Section 
6  below,  this  can  be  explained  by  the  size  of  the  throats  selected 
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of  liquid  saturation  for  various  fractions  of  hydrophilic  pores  (solid  line  for  a 


during  the  invasion  of  the  network  by  the  liquid.  Narrower  throats 
are  favoured  in  a  hydrophilic  network  whereas  the  invasion  occurs 
preferentially  in  wider  throats  in  a  hydrophobic  network. 

Experimental  measurements  of  relative  permeabilities  are 
reported  in  [29]  for  a  GDL  treated  with  a  5  wt%  PTFE  solution  and 
in  [30],  but  only  for  untreated  (non-teflonized)  materials.  Our  sim¬ 
ulated  values  of  the  gas  phase  relative  permeability  are  somewhat 
greater  than  the  experimental  data  reported  in  [29]  but  in  reason¬ 
able  agreement  with  the  lattice  Boltzmann  computation  of  I<rg  also 
reported  in  [29].  More  accurate  experimental  data  would  be  nec¬ 
essary  to  make  a  more  informative  comparison.  Our  results  for  a 
sufficiently  hydrophobic  network  are  also  in  good  agreement  with 
the  lattice  Boltzmann  computation  of  Krg  and  I<rt  reported  in  [14] 
for  a  small  volume  of  a  model  hydrophobic  GDL. 

The  measured  values  of  through-plane  gas  phase  relative  per¬ 
meability  reported  in  [30]  are  found  to  be  very  low,  much  lower 
than  predicted  with  the  PNM  or  reported  in  [14,29].  According 
to  [30]  this  low  values  are  attributed  to  an  experimental  artefact. 
However,  the  GDLs  being  untreated  in  this  experiment,  compact 
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invasion  patterns  are  possible,  see  Fig.  2  and  Ref.  [3],  and  might 
also  explain  in  part  the  results  reported  in  [30]. 

6.  Effective  binary  diffusion  coefficient 

The  (through-plane)  effective  binary  diffusion  coefficient  De^  is 
an  interesting  parameter  because  it  reflects  the  effect  of  the  pore 
blockage  by  the  water  on  the  gas  transport  through  the  GDL.  In  fuel 
cell  modelling  literature,  Dejy(S)  is  often  represented  by  expression 
of  the  form,  e.g.  [17] 

Deff{S)  =  e1'5(l  -  S)1,5D  (5) 

where  D  is  the  binary  molecular  diffusion  coefficient.  Here  we  use 
the  pore  network  model  to  compute  Deff  and  it  is  shown  that  Eq.  (5) 
does  not  lead  to  a  good  match  with  the  pore  network  simulations. 

The  procedure  for  computing  the  effective  binary  diffusion  coef¬ 
ficient  from  pore  network  simulations  is  analogous  to  the  one  for 
the  effective  permeability  and  is  described  for  example  in  [17].  For 
each  stage  of  pore  network  occupancy  by  the  water  obtained  as 
described  in  Section  4,  local  diffusive  conductances  are  assigned  to 
each  pore  and  throat  occupied  by  the  gas  phase  and  a  linear  sys¬ 
tem  of  equations  is  formed  by  expressing  the  mass  conservation 
at  each  gaseous  pore  of  the  network  and  imposing  a  concentra¬ 
tion  difference  across  the  network.  Solving  numerically  the  linear 
system  give  the  concentration  of  the  considered  species  (i.e.  oxy¬ 
gen)  at  each  gaseous  pore  of  network.  This  allows  to  compute  the 
diffusive  flux  across  the  network  and  then  to  determine  the  effec¬ 
tive  diffusion  coefficient  from  the  mathematical  expression  of  the 
through-plane  macroscopic  flux  J ,  which  is  expressed  as 

j  An  ( Qnlet  —  Qiutlet 

J  =  ADeff  ^ - - - 

where  A  the  cross-section  area  of  the  porous  medium,  L  is  the 
porous  medium  thickness  and  Cinlet  (resp.  Coutlet)  is  the  concen¬ 
tration  imposed  at  the  inlet  (resp.  outlet).  The  procedure  used  here 
is  exactly  the  same  as  the  one  described  in  [17]  and  therefore  the 
details  are  not  given  again  for  the  sake  of  brevity. 

The  simulation  results  are  presented  in  Fig.  5.  They  are  not  sur¬ 
prisingly  in  good  agreement  with  the  PNM  simulations  presented 
in  [17]  for  a  fully  hydrophobic  network  since  the  computational 
procedure  is  the  same  as  in  [17]  (for  the  case  referred  to  as  “case  1” 
in  [1 7]).  More  interesting  is  the  impact  of  the  fraction  of  hydrophilic 
pores  /.  For  a  given  intermediate  saturation,  Dejj  increases  signifi¬ 
cantly  with/ for />/c,  i.e./>50%  in  Fig.  5.  Thus  a  rapid  conclusion 
would  be  that  a  sufficiently  hydrophilic  GDL  is  better  since  the 
effective  diffusion  coefficient  increases  with  /  for  f>fc.  However, 
one  should  keep  in  mind  that  the  phase  distribution  evolution  is 
obtained  here  as  described  in  Section  4,  i.e.  as  the  result  of  suc¬ 
cessive  states  of  hydrostatic  equilibrium  corresponding  to  small 
increment  in  the  pressure  of  the  water  dPw.  This  procedure  leads 
to  the  preferential  invasion  of  the  smallest  throats  and  pores  in 
a  hydrophilic  network  and  to  the  preferential  invasion  of  largest 
pores  in  a  hydrophobic  networks  (this  can  be  easily  seen  from  the 
expressions  of  capillary  pressure  thresholds  in  Section  3).  There¬ 
fore,  for  a  given  saturation,  the  fraction  of  the  pore  space  left  free 
for  the  gas  transport  is  primarily  made  of  large  pores  and  throats  in  a 
hydrophilic  network  and  of  small  pores  and  throats  in  a  hydropho¬ 
bic  network.  This  explains  why  the  effective  diffusion  coefficient 
is  greater  in  the  hydrophilic  network  and  increases  with/ for  f>fc. 
The  same  explanation  holds  as  regards  the  evolution  of  the  relative 
permeabilities  reported  in  Section  5. 

As  presented  in  the  next  sections,  the  conclusion  is  exactly  the 
contrary,  i.e.  a  too  hydrophilic  GDL  is  detrimental  to  gas  trans¬ 
port,  when  water  invades  the  GDL  according  to  the  invasion  rules 
described  in  Section  7.  As  explained  in  Section  4,  the  procedure  used 
in  Sections  4-6  to  compute  the  invasion  of  the  network  by  the  liquid 
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Fig.  5.  Evolution  of  through-plane  binary  effective  diffusion  coefficient  as  a  func¬ 
tion  of  liquid  saturation  for  various  fractions  of  hydrophilic  pores  (solid  line  for  a 
hydrophilic  pore  fraction  below  the  percolation  threshold  fc,  dashed  line  above /c): 
(a)  fractions  of  hydrophilic  pores  below  the  percolation  threshold,  (b)  fractions  of 
hydrophilic  pores  above  the  percolation  threshold. 

water  corresponds  to  the  typical  procedure  used  to  measure  a  cap¬ 
illary  pressure  curve.  In  this  case,  the  water  invasion  is  controlled 
by  increasing  step  by  step  the  pressure  difference  between  the  two 
fluids.  The  completely  different  invasion  condition  used  in  the  next 
sections  corresponds  to  imposing  a  small  invasion  flow  rate.  This 
condition  is  supposedly  much  more  representative  of  water  inva¬ 
sion  in  a  GDL  in  a  fuel  cell  than  the  pressure  difference  controlled 
invasion.  This  illustrates  that  care  should  be  exercised  when  ana¬ 
lyzing  two-phase  flows  in  this  type  of  system  or  when  devising  ex 
situ  experiments  for  characterizing  GDL  transport  properties. 

As  can  be  seen  from  Fig.  5,  Eq.  (5)  overpredicts  Dejf  (5)  sig¬ 
nificantly.  Thus  using  Eq.  (5)  leads  to  underestimating  the  mass 
transport  limitation  in  the  gas  diffusion  layer.  The  comparisons 
with  the  experimental  results  reported  in  [31]  (for  a  dry  GDL) 
indicate  that  this  holds  for  all  the  existing  effective  diffusion  theo¬ 
retical  models.  Compared  to  the  experimental  results  presented  in 
[31  ]  (Dejff(0)/D  ~  0.32  for  an  untreated  GDL),  the  PNM  also  overpre¬ 
dicts  Deff  with  DeJff(0)/D  ~  0.5,  but  less  than  the  available  theoretical 
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models  (Dejy{0)/D~0.6-0.7,  see  [31]).  Although  the  experimental 
evaluation  of  Deff  is  not  easy  due  to  the  thin  nature  of  a  GDL  and 
therefore  additional  measurements  are  desirable,  this  can  be  again 
an  indication  that  better  pore  network  models  of  fibrous  materials 
need  to  be  developed.  In  passing,  we  note  that  the  experimental 
results  reported  in  [31]  indicate  a  significant  influence  of  Teflon 
treatment  on  Defj{ 0).  This  would  indicate  that  the  Teflon  treatment 
modifies  the  pore  space  sufficiently  for  affecting  the  effective  diffu¬ 
sion  coefficient  but  not,  according  to  [1  ],  for  affecting  the  capillary 
pressure  curve.  This  is  somewhat  surprising  and  would  need  fur¬ 
ther  clarification.  As  mentioned  before,  this  effect  of  Teflon  on  pore 
size  and  pore  volume  is  not  taken  into  account  in  our  simulations. 

7.  Quasi-static  invasion 

In  the  preceding  sections,  the  PNM  was  essentially  used  to  com¬ 
pute  some  parameters  of  the  classical  phenomenological  model 
of  two-phase  flow  in  porous  media.  This  is  a  somewhat  obvious 
objective  since  most  of  the  simulations  of  two-phase  flow  in  fuel 
cells  are  based  on  this  classical  model.  As  discussed  for  example  in 
[22],  there  are,  however,  serious  concerns  about  using  the  classi¬ 
cal  phenomenological  model  to  simulate  water  invasion  in  a  thin 
system  like  a  GDL.  There  is  a  lack  of  length  separation  (a  GDL  is 
typically  a  few  pore  sizes  thick)  and  the  dominant  two  phase  flow 
regime  in  a  hydrophobic  GDL  is  expected  to  be  the  capillary  fin¬ 
gering  regime,  a  regime  which  leads  to  fractal  distributions  of  the 
fluids  within  the  pore  space  and  which  is  therefore  not  compatible 
with  the  traditional  concept  of  R.E.V.  (representative  elementary 
volume)  underlying  the  classical  continuum  approach.  As  a  result, 
the  classical  continuum  model  is  likely  to  lead  to  poor  approxima¬ 
tions  of  the  real  fields  and  this  is  illustrated  in  [22].  Because  of  the 
weaknesses  of  the  classical  model,  it  can  be  much  more  instructive 
to  use  the  PNM  for  predicting  directly  the  water  invasion  pattern 
within  the  GDL. 

As  mentioned  before,  water  invasion  in  a  GDL  for  the  operat¬ 
ing  conditions  of  a  PEMFC  is  expected  to  be  dominated  by  capillary 
effects  and  therefore  can  be  simulated  as  a  quasi-static  invasion 
process.  As  briefly  evoked  in  the  preceding  section,  Section  6,  one 
can  distinguish  two  main  types  of  boundary  condition  (B.C.)  to  sim¬ 
ulate  a  quasi-static  liquid  water  invasion:  (i)  pressure  B.C.  and  (ii) 
flow  rate  B.C.  As  described  in  Section  4,  the  pressure  B.C.  consists 
in  increasing  step  by  step  the  pressure  difference  between  liquid 
water  and  the  gas  phase.  The  flow  rate  B.C.  consists  in  assuming 
that  liquid  water  invades  the  system  as  a  result  of  a  small  imposed 
flow  rate  at  the  inlet.  In  this  case,  it  is  classically  assumed  that 
all  the  throats  at  the  inlet  of  the  network  are  connected  to  a  con¬ 
tinuous  liquid  reservoir  at  the  inlet  and  that  liquid  is  very  slowly 
injected  at  a  specified  flow  rate  in  this  reservoir.  This  increases  very 
gently  the  pressure  (supposed  uniform)  in  the  reservoir  until  the 
invasion  of  the  adjacent  porous  medium  can  start.  Then  the  pres¬ 
sure  in  the  reservoir  typically  fluctuates  around  a  mean  value  as 
the  water  invasion  continues  in  the  porous  medium  up  to  break¬ 
through.  Starting  from  a  dry  network,  the  corresponding  invasion 
algorithm  reads, 

(1)  Identify  the  element  (pore  or  throat)  adjacent  to  the  already 
invaded  region  that  has  the  lowest  invasion  capillary  pressure 
threshold  (see  Section  3). 

(2)  Invade  the  identified  element  in  1  and  update  the  phase  distri¬ 
bution  in  the  network.  Return  to  (1 ). 

The  invasion  stops  when  water  reaches  the  GDL  outlet.  This 
stage  of  the  invasion  is  referred  to  as  the  breakthrough  (BT). 

As  discussed  in  [22]  or  [32],  there  is  still,  however,  a  question 
about  the  flow  rate  boundary  condition  to  impose  at  the  inlet  of  the 


f(%) 


Fig.  6.  (a)  Evolution  of  overall  liquid  saturation  at  breakthrough  as  a  function  of  the 
fraction  of  hydrophilic  pores;  (b)  percolation  probability  as  a  function  of  the  fraction 
of  hydrophilic  pores. 

GDL.  As  discussed  above,  the  classical  flow  rate  boundary  condition 
is  to  assume  that  all  throats  of  the  network  at  the  inlet  are  connected 
to  a  water  reservoir.  This  condition  can  be  representative  of  flow 
rate  controlled  ex  situ  experiments,  e.g.  [33  ].  For  a  GDL  in  a  fuel  cell, 
the  boundary  condition  can  be  expected  to  be  different  since  water 
comes  not  from  a  reservoir  but  from  adjacent  porous  layers  (MPL 
and  catalyst  layer).  As  discussed  in  [32],  it  is  envisioned  that  water 
enters  rather  the  GDL  through  a  series  of  independent  injection 
points.  As  shown  in  [32],  the  boundary  condition  has  a  significant 
impact  on  the  GDL  water  occupancy.  Thus  it  would  be  interesting  to 
consider  the  multiple  injection  points  boundary  condition  as  well. 
For  the  sake  of  brevity,  we  only  consider,  however,  the  reservoir 
type  boundary  condition.  The  effect  of  the  multiple  injection  points 
boundary  condition  will  be  explored  in  a  future  work. 

The  evolution  of  the  overall  liquid  water  saturation  S  at  BT  as  a 
function  of  the  fraction  of  hydrophilic  pores  /  is  shown  in  Fig.  6a. 
As  in  any  classical  percolation  problem  [34],  the  hydrophilic  pores 
and  throats  form  isolated  hydrophilic  pore  and  throat  clusters  for 
sufficiently  low  values  off.  When/is  sufficiently  large,  there  exists 
a  cluster  of  hydrophilic  pores  and  throats  spanning  the  system.  The 
critical  hydrophilic  pore  fraction  below  which  such  a  cluster  does 
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Fig.  7.  Liquid  saturation  profile  along  the  thickness  at  breakthrough  for  various 
fractions  of  hydrophilic  pores  (solid  line  for  a  hydrophilic  pore  fraction  below  the 
percolation  threshold  fc,  dashed  line  above  fc). 

not  exist  is  the  percolation  threshold  /.  In  a  finite  size  system,  the 
percolation  transition  is  not  sharp  and  some  realizations  are  perco¬ 
lating  below  /  whereas  others  are  percolating  for  values  of/greater 
than  /.  This  is  illustrated  in  Fig.  6b,  which  shows  the  evolution  of 
the  percolation  probability  F  as  a  function  of/ (for  a  given  value 
/,  F  is  the  fraction  of  percolating  realizations  among  the  100  gen¬ 
erated  realizations).  The  percolation  threshold  corresponds  to  the 
rapid  transition  in  the  evolution  of  Fin  Fig.  6b.  As  can  be  seen  from 
Fig.  6b,  fc  is  of  the  order  of  46-48%  for  our  system. 

Since  liquid  water  favours  hydrophilic  elements,  liquid  water 
invasion  takes  place  through  a  path  of  hydrophilic  pores  and  throats 
above  /  and  the  hydrophobic  pores  and  throats  are  left  free  of 
water,  i.e.  available  for  gas  transport.  Below  /,  some  hydropho¬ 
bic  pores  and  throats  are  necessarily  invaded  since  the  hydrophilic 
pores  and  throats  do  not  form  a  percolating  cluster.  However, 
a  hydrophilic  cluster  is  systematically  fully  invaded  when  liquid 
water  meets  such  a  hydrophilic  cluster  during  the  invasion.  We 
know  from  Section  3  and  Ref.  [3]  that  the  invasion  pattern  and  the 
local  invasion  mechanisms  are  significantly  different  for  0  «  80°  and 
0^110°  (as  illustrated  also  in  Fig.  2).  Interestingly,  this  does  not 
have,  however,  a  great  impact  on  the  overall  saturation  at  break¬ 
through  over  a  large  range  of/,  i.e.  as  long  as /</  as  shown  in  Fig.  6a. 
This  can  be  qualitatively  understood  as  follows.  For/close  to  0,  the 
water  invasion  process  is  essentially  an  invasion  percolation  (IP) 
process  as  discussed  in  Section  3.  As  it  is  well  known,  the  subset  of 
pores  invaded  at  breakthrough  as  a  result  of  the  IP  process  is  very 
similar  to  an  ordinary  percolation  cluster,  i.e.  the  percolating  cluster 
formed  by  the  hydrophilic  pores  and  throats  at/c.  Hence,  since  the 
pore  size  distribution  is  generally  narrow  in  a  GDL,  we  expect  that 
S(0)&S(fc).  For  intermediate  /between  0  and  /,  the  invaded  pore 
subset  can  be  envisioned  as  a  chain  of  invaded  hydrophilic  pore 
clusters  connected  by  paths  of  hydrophobic  pores,  i.e.  a  structure 
close  to  a  percolation  cluster  at  /.  The  saturation  increases  rapidly 
above/c.  The  minimum  saturation  in  Fig.  6a  is  observed  for  a  fully 
hydrophobic  network  (/=  0).  Interestingly,  a  local  minimum  satu¬ 
ration  is  observed  in  Fig.  6a  for  /«/.  The  reasons  for  the  existence 
of  this  local  minimum  are,  however,  unclear  and  this  is  left  for  a 
future  work. 

The  main  practical  conclusion  from  the  results  shown  in  Fig.  6, 
is  that  the  PTFE  load  (assuming  here  that  there  is  a  direct  link 
between  the  PTFE  load  and  the  fraction  of  hydrophobic  pores)  has 


Fig.  8.  Evolution  of  through-plane  binary  effective  diffusion  coefficient  at  break¬ 
through  as  a  function  of  the  fraction  of  hydrophilic  pores. 

little  impact  on  the  water  management  problem  as  long  as  the  frac¬ 
tion  of  hydrophilic  pores  and  throats  is  not  significantly  larger  than 
the  percolation  threshold  of  the  system. 

This  is  further  illustrated  with  the  saturation  profiles  at  BT 
depicted  in  Fig.  7.  These  profiles  are  obtained  by  computing  the 
saturation  over  successive  slices  of  network  in  the  through-plane 
direction  (10  slices  for  our  40  x  40  x  10  network)  for  each  real¬ 
ization.  Again  the  profile  is  not  very  sensitive  to  the  fraction  of 
hydrophilic  pores  as  long  as  /</  except  near  the  inlet  where  the 
minimum  saturation  in  the  first  slice  (S  ~  0.35)  is  obtained  for  a  fully 
hydrophobic  network  whereas  S~0.43  for  /approaching/  (curve 
for/=  40%  in  Fig.  7).  As  can  be  seen  from  Fig.  7,  the  pore  blockage  by 
the  water  at  BT  increases  significantly  for />/.  Again,  these  results 
suggest  that  a  loss  in  hydrophobicity  in  our  model  system  (assum¬ 
ing  for  instance  a  fully  hydrophobic  network  initially)  should  not 
lead  to  a  discernable  decrease  in  the  fuel  cell  performance  in  a  first 
phase  (i.e.  as  long  as /increases  below/)  and  then  to  a  significant 
loss  in  performance  when  the  GDL  becomes  sufficiently  hydrophilic 
(when/increases  above/).  The  evolution  of  the  effective  diffusion 
coefficient  at  breakthrough  presented  in  the  next  section  confirms 
this  conclusion. 

8.  Effective  diffusion  coefficient  and  relative  permeabilities 
at  breakthrough 

The  evolution  of  the  effective  diffusion  coefficient  (through- 
plane  coefficient  over  the  whole  thickness  of  the  network)  at  BT  as 
a  function  of  hydrophilic  pore  fraction  is  shown  in  Fig.  8  together 
with  the  evolution  of  the  overall  saturation.  The  evolution  given 
by  Eq.  (5)  is  also  plotted  in  Fig.  8  and  can  be  seen  as  overpredict¬ 
ing  our  results  (see  the  end  of  Section  6  in  this  respect).  Again,  one 
can  first  distinguish  two  main  domains.  For/below/c,  Dejj  remains 
high  (of  the  order  of  0.4-0.5D)  whereas  it  decreases  steadily  up 
to  0  (which  is  reached  for /~ 85%)  as /increases  above/.  As  can 
be  seen  from  Fig.  8,  the  maximum  value  of  Deff  is  obtained  for  a 
fully  hydrophobic  network,  i.e./=0.  A  local  maximum  in  the  evo¬ 
lution  of  Deff  is  observed  for  /~/c,  right  before  the  beginning  of  the 
steady  drop  above  /.  In  summary,  these  results  suggest  that  the 
loss  in  performance  due  a  diminished  access  of  gas  to  the  catalyst 
layer  through  the  GDL  resulting  from  a  change  in  wettability  is  only 
marked  when />/.  Below/,  the  best  performance  is  expected  for 
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Fig.  9.  Evolution  of  through-plane  relative  permeabilities  at  breakthrough  as  a  func¬ 
tion  of  the  fraction  of  hydrophilic  pores. 

a  fully  hydrophobic  GDL  with  also  good  performances  expected  for 
/in  the  range  [0-/c].  Of  course,  these  conclusions  are  for  our  model 
GDL  The  fact  that  a  GDL  becomes  or  not  a  mixed  wettability  system 
as  a  result  of  aging  or  due  to  the  temperature  conditions  prevail¬ 
ing  in  an  operating  PEMFC  remains  an  open  question.  Also,  even  if 
the  mixed  wettability  assumption  is  correct,  there  is  no  evidence  at 
the  moment  that  the  GDL  can  be  considered  as  a  system  of  purely 
random  wettability  as  considered  in  this  paper.  Nevertheless,  our 
simple  model  is  consistent  with  the  expected  effect  of  a  sufficient 
loss  in  hydrophobicity  and  also  suggests  that  a  mixed  wettability 
GDL  can  be  an  efficient  solution  as  regards  the  water  management 
problem  provided  that  the  hydrophobic  fraction  of  the  pore  space  is 
sufficiently  important.  We  note  that  a  somewhat  similar  conclusion 
was  reached  in  [21]. 

Similar  conclusion  can  be  obtained  from  the  evolution  of  Kr  at 
breakthrough  shown  in  Fig.  9. 

9.  Conclusions 

In  this  paper,  we  explored  the  effect  of  a  change  in  wettability 
from  pore  network  simulations  for  a  model  GDL  considered  as  a 
cubic  pore  network  of  random  mixed  wettability.  The  study  illus¬ 
trates  well  the  possibilities  offered  by  pore  network  simulations  for 
establishing  the  link  between  a  change  in  the  local  properties  and 
a  change  in  the  macroscopic  equilibrium  (capillary  pressure  curve) 
or  transport  properties  of  the  porous  medium. 

The  simulations  lead  to  results  compatible  with  the  most  recent 
measurements  of  the  water  invasion  capillary  pressure  curve  for 
GDLs  treated  with  PTFE.  A  somewhat  surprising  result  is  that  the 
macroscopic  properties  studied  are  independent  of  the  fraction  of 
hydrophilic  pores  and  throats  in  the  system  as  long  as  this  fraction 
is  below  a  certain  threshold  corresponding  in  our  model  to  the  per¬ 
colation  threshold  of  the  hydrophilic  pore  network.  This  result  is 
somewhat  analogous  to  the  experimental  results  showing  that  the 
capillary  pressure  curve  is  essentially  independent  of  the  PTFE  load¬ 
ing  as  long  as  the  loading  is  equal  to  or  greater  than  5  wt%.  Thus  it  is 
difficult  to  infer  from  the  available  experimental  results  concerning 
the  capillary  pressure  curve  of  GDLs  whether  it  is  more  appropriate 
to  consider  a  GDL  as  a  system  of  spatially  (quasi-)  uniform  wetta¬ 
bility  (as  suggested  from  the  analysis  and  results  reported  in  [1] 
for  example)  or  rather  as  a  system  of  mixed  wettability  (i.e.  locally 
hydrophilic  in  some  regions  of  the  pore  space  and  hydrophobic 
elsewhere). 


The  same  type  of  question,  i.e.  a  spatially  uniform  change  or  on 
the  contrary  local  changes,  is  still  open  as  regards  the  change  in 
wettability  that  might  be  due  to  aging  or  temperature  effect.  We 
explored  the  option  of  a  change  in  wettability  due  to  an  increase  in 
the  fraction  of  hydrophilic  pores  and  throats,  assuming  moreover 
a  random  distribution  of  the  hydrophilic  pores. 

This  model  leads  to  delineate  two  main  domains,  below  and 
above  the  hydrophilic  pore  percolation  threshold.  Below  the  per¬ 
colation  threshold,  all  the  investigated  properties  of  the  system  are 
essentially  insensitive  to  a  change  in  the  fraction  of  hydrophilic 
pores.  Above  the  percolation  threshold,  the  effect  of  a  change  in 
the  fraction  of  hydrophilic  pores  becomes  quite  significant  and  the 
gas  access  through  the  GDL  becomes  quite  poor  when  the  fraction 
of  hydrophilic  pores  is  high.  The  identification  of  these  two  regimes 
may  help  interpret  experimental  evidences  of  fuel  cell  performance 
loss  due  to  change  in  wettability. 

We  note,  however,  that  the  other  “extreme”  scenario,  that  is  a 
progressive  and  spatially  uniform  change  in  the  wettability  (leading 
to  water  distributions  from  capillary  fingering  patterns  to  com¬ 
pact  ones)  is  difficult  to  study  from  pore  network  simulations  at 
the  moment  because  there  is  no  three-dimensional  version  of  PNM 
allowing  to  study  accurately  the  evolution  of  invasion  patterns  in 
the  contact  angle  intermediate  range  [70-120°]  of  interest  for  a 
GDL.  So  we  do  not  know  if  the  second  scenario  of  changing  uni¬ 
form  wettability  would  lead  to  an  evolution  of  the  macroscopic 
properties  markedly  different  from  the  mixed  wettability  scenario 
explored  in  this  paper.  Hence,  further  work  on  two-phase  flow 
in  the  intermediate  range  of  contact  angle  characterizing  GDLs  is 
needed. 

Finally,  we  note  that  the  simulations  presented  in  this  paper 
suggest  that  the  best  performances  can  be  expected  not  only  with 
a  fully  hydrophobic  GDL  but  also  with  partially  hydrophilic  GDLs. 
Hence  a  partial  teflonization  of  the  GDL  is  not  necessarily  a  problem 
provided  that  the  GDL  contains  a  sufficient  fraction  of  hydrophobic 
pores. 
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